home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / blas / zher2k.f < prev    next >
Text File  |  1997-06-25  |  14KB  |  373 lines
  1.       SUBROUTINE ZHER2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, BETA,
  2.      $                   C, LDC )
  3. *     .. Scalar Arguments ..
  4.       CHARACTER          TRANS, UPLO
  5.       INTEGER            K, LDA, LDB, LDC, N
  6.       DOUBLE PRECISION   BETA
  7.       COMPLEX*16         ALPHA
  8. *     ..
  9. *     .. Array Arguments ..
  10.       COMPLEX*16         A( LDA, * ), B( LDB, * ), C( LDC, * )
  11. *     ..
  12. *
  13. *  Purpose
  14. *  =======
  15. *
  16. *  ZHER2K  performs one of the hermitian rank 2k operations
  17. *
  18. *     C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) + beta*C,
  19. *
  20. *  or
  21. *
  22. *     C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A + beta*C,
  23. *
  24. *  where  alpha and beta  are scalars with  beta  real,  C is an  n by n
  25. *  hermitian matrix and  A and B  are  n by k matrices in the first case
  26. *  and  k by n  matrices in the second case.
  27. *
  28. *  Parameters
  29. *  ==========
  30. *
  31. *  UPLO   - CHARACTER*1.
  32. *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
  33. *           triangular  part  of the  array  C  is to be  referenced  as
  34. *           follows:
  35. *
  36. *              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
  37. *                                  is to be referenced.
  38. *
  39. *              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
  40. *                                  is to be referenced.
  41. *
  42. *           Unchanged on exit.
  43. *
  44. *  TRANS  - CHARACTER*1.
  45. *           On entry,  TRANS  specifies the operation to be performed as
  46. *           follows:
  47. *
  48. *              TRANS = 'N' or 'n'    C := alpha*A*conjg( B' )          +
  49. *                                         conjg( alpha )*B*conjg( A' ) +
  50. *                                         beta*C.
  51. *
  52. *              TRANS = 'C' or 'c'    C := alpha*conjg( A' )*B          +
  53. *                                         conjg( alpha )*conjg( B' )*A +
  54. *                                         beta*C.
  55. *
  56. *           Unchanged on exit.
  57. *
  58. *  N      - INTEGER.
  59. *           On entry,  N specifies the order of the matrix C.  N must be
  60. *           at least zero.
  61. *           Unchanged on exit.
  62. *
  63. *  K      - INTEGER.
  64. *           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
  65. *           of  columns  of the  matrices  A and B,  and on  entry  with
  66. *           TRANS = 'C' or 'c',  K  specifies  the number of rows of the
  67. *           matrices  A and B.  K must be at least zero.
  68. *           Unchanged on exit.
  69. *
  70. *  ALPHA  - COMPLEX*16         .
  71. *           On entry, ALPHA specifies the scalar alpha.
  72. *           Unchanged on exit.
  73. *
  74. *  A      - COMPLEX*16       array of DIMENSION ( LDA, ka ), where ka is
  75. *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
  76. *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
  77. *           part of the array  A  must contain the matrix  A,  otherwise
  78. *           the leading  k by n  part of the array  A  must contain  the
  79. *           matrix A.
  80. *           Unchanged on exit.
  81. *
  82. *  LDA    - INTEGER.
  83. *           On entry, LDA specifies the first dimension of A as declared
  84. *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
  85. *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
  86. *           be at least  max( 1, k ).
  87. *           Unchanged on exit.
  88. *
  89. *  B      - COMPLEX*16       array of DIMENSION ( LDB, kb ), where kb is
  90. *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
  91. *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
  92. *           part of the array  B  must contain the matrix  B,  otherwise
  93. *           the leading  k by n  part of the array  B  must contain  the
  94. *           matrix B.
  95. *           Unchanged on exit.
  96. *
  97. *  LDB    - INTEGER.
  98. *           On entry, LDB specifies the first dimension of B as declared
  99. *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
  100. *           then  LDB must be at least  max( 1, n ), otherwise  LDB must
  101. *           be at least  max( 1, k ).
  102. *           Unchanged on exit.
  103. *
  104. *  BETA   - DOUBLE PRECISION            .
  105. *           On entry, BETA specifies the scalar beta.
  106. *           Unchanged on exit.
  107. *
  108. *  C      - COMPLEX*16          array of DIMENSION ( LDC, n ).
  109. *           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
  110. *           upper triangular part of the array C must contain the upper
  111. *           triangular part  of the  hermitian matrix  and the strictly
  112. *           lower triangular part of C is not referenced.  On exit, the
  113. *           upper triangular part of the array  C is overwritten by the
  114. *           upper triangular part of the updated matrix.
  115. *           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
  116. *           lower triangular part of the array C must contain the lower
  117. *           triangular part  of the  hermitian matrix  and the strictly
  118. *           upper triangular part of C is not referenced.  On exit, the
  119. *           lower triangular part of the array  C is overwritten by the
  120. *           lower triangular part of the updated matrix.
  121. *           Note that the imaginary parts of the diagonal elements need
  122. *           not be set,  they are assumed to be zero,  and on exit they
  123. *           are set to zero.
  124. *
  125. *  LDC    - INTEGER.
  126. *           On entry, LDC specifies the first dimension of C as declared
  127. *           in  the  calling  (sub)  program.   LDC  must  be  at  least
  128. *           max( 1, n ).
  129. *           Unchanged on exit.
  130. *
  131. *
  132. *  Level 3 Blas routine.
  133. *
  134. *  -- Written on 8-February-1989.
  135. *     Jack Dongarra, Argonne National Laboratory.
  136. *     Iain Duff, AERE Harwell.
  137. *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  138. *     Sven Hammarling, Numerical Algorithms Group Ltd.
  139. *
  140. *  -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1.
  141. *     Ed Anderson, Cray Research Inc.
  142. *
  143. *
  144. *     .. External Functions ..
  145.       LOGICAL            LSAME
  146.       EXTERNAL           LSAME
  147. *     ..
  148. *     .. External Subroutines ..
  149.       EXTERNAL           XERBLA
  150. *     ..
  151. *     .. Intrinsic Functions ..
  152.       INTRINSIC          DBLE, DCONJG, MAX
  153. *     ..
  154. *     .. Local Scalars ..
  155.       LOGICAL            UPPER
  156.       INTEGER            I, INFO, J, L, NROWA
  157.       COMPLEX*16         TEMP1, TEMP2
  158. *     ..
  159. *     .. Parameters ..
  160.       DOUBLE PRECISION   ONE
  161.       PARAMETER          ( ONE = 1.0D+0 )
  162.       COMPLEX*16         ZERO
  163.       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  164. *     ..
  165. *     .. Executable Statements ..
  166. *
  167. *     Test the input parameters.
  168. *
  169.       IF( LSAME( TRANS, 'N' ) ) THEN
  170.          NROWA = N
  171.       ELSE
  172.          NROWA = K
  173.       END IF
  174.       UPPER = LSAME( UPLO, 'U' )
  175. *
  176.       INFO = 0
  177.       IF( ( .NOT.UPPER ) .AND. ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN
  178.          INFO = 1
  179.       ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ) .AND.
  180.      $         ( .NOT.LSAME( TRANS, 'C' ) ) ) THEN
  181.          INFO = 2
  182.       ELSE IF( N.LT.0 ) THEN
  183.          INFO = 3
  184.       ELSE IF( K.LT.0 ) THEN
  185.          INFO = 4
  186.       ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
  187.          INFO = 7
  188.       ELSE IF( LDB.LT.MAX( 1, NROWA ) ) THEN
  189.          INFO = 9
  190.       ELSE IF( LDC.LT.MAX( 1, N ) ) THEN
  191.          INFO = 12
  192.       END IF
  193.       IF( INFO.NE.0 ) THEN
  194.          CALL XERBLA( 'ZHER2K', INFO )
  195.          RETURN
  196.       END IF
  197. *
  198. *     Quick return if possible.
  199. *
  200.       IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
  201.      $    ( BETA.EQ.ONE ) ) )RETURN
  202. *
  203. *     And when  alpha.eq.zero.
  204. *
  205.       IF( ALPHA.EQ.ZERO ) THEN
  206.          IF( UPPER ) THEN
  207.             IF( BETA.EQ.DBLE( ZERO ) ) THEN
  208.                DO 20 J = 1, N
  209.                   DO 10 I = 1, J
  210.                      C( I, J ) = ZERO
  211.    10             CONTINUE
  212.    20          CONTINUE
  213.             ELSE
  214.                DO 40 J = 1, N
  215.                   DO 30 I = 1, J - 1
  216.                      C( I, J ) = BETA*C( I, J )
  217.    30             CONTINUE
  218.                   C( J, J ) = BETA*DBLE( C( J, J ) )
  219.    40          CONTINUE
  220.             END IF
  221.          ELSE
  222.             IF( BETA.EQ.DBLE( ZERO ) ) THEN
  223.                DO 60 J = 1, N
  224.                   DO 50 I = J, N
  225.                      C( I, J ) = ZERO
  226.    50             CONTINUE
  227.    60          CONTINUE
  228.             ELSE
  229.                DO 80 J = 1, N
  230.                   C( J, J ) = BETA*DBLE( C( J, J ) )
  231.                   DO 70 I = J + 1, N
  232.                      C( I, J ) = BETA*C( I, J )
  233.    70             CONTINUE
  234.    80          CONTINUE
  235.             END IF
  236.          END IF
  237.          RETURN
  238.       END IF
  239. *
  240. *     Start the operations.
  241. *
  242.       IF( LSAME( TRANS, 'N' ) ) THEN
  243. *
  244. *        Form  C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) +
  245. *                   C.
  246. *
  247.          IF( UPPER ) THEN
  248.             DO 130 J = 1, N
  249.                IF( BETA.EQ.DBLE( ZERO ) ) THEN
  250.                   DO 90 I = 1, J
  251.                      C( I, J ) = ZERO
  252.    90             CONTINUE
  253.                ELSE IF( BETA.NE.ONE ) THEN
  254.                   DO 100 I = 1, J - 1
  255.                      C( I, J ) = BETA*C( I, J )
  256.   100             CONTINUE
  257.                   C( J, J ) = BETA*DBLE( C( J, J ) )
  258.                ELSE
  259.                   C( J, J ) = DBLE( C( J, J ) )
  260.                END IF
  261.                DO 120 L = 1, K
  262.                   IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) )
  263.      $                 THEN
  264.                      TEMP1 = ALPHA*DCONJG( B( J, L ) )
  265.                      TEMP2 = DCONJG( ALPHA*A( J, L ) )
  266.                      DO 110 I = 1, J - 1
  267.                         C( I, J ) = C( I, J ) + A( I, L )*TEMP1 +
  268.      $                              B( I, L )*TEMP2
  269.   110                CONTINUE
  270.                      C( J, J ) = DBLE( C( J, J ) ) +
  271.      $                           DBLE( A( J, L )*TEMP1+B( J, L )*TEMP2 )
  272.                   END IF
  273.   120          CONTINUE
  274.   130       CONTINUE
  275.          ELSE
  276.             DO 180 J = 1, N
  277.                IF( BETA.EQ.DBLE( ZERO ) ) THEN
  278.                   DO 140 I = J, N
  279.                      C( I, J ) = ZERO
  280.   140             CONTINUE
  281.                ELSE IF( BETA.NE.ONE ) THEN
  282.                   DO 150 I = J + 1, N
  283.                      C( I, J ) = BETA*C( I, J )
  284.   150             CONTINUE
  285.                   C( J, J ) = BETA*DBLE( C( J, J ) )
  286.                ELSE
  287.                   C( J, J ) = DBLE( C( J, J ) )
  288.                END IF
  289.                DO 170 L = 1, K
  290.                   IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) )
  291.      $                 THEN
  292.                      TEMP1 = ALPHA*DCONJG( B( J, L ) )
  293.                      TEMP2 = DCONJG( ALPHA*A( J, L ) )
  294.                      DO 160 I = J + 1, N
  295.                         C( I, J ) = C( I, J ) + A( I, L )*TEMP1 +
  296.      $                              B( I, L )*TEMP2
  297.   160                CONTINUE
  298.                      C( J, J ) = DBLE( C( J, J ) ) +
  299.      $                           DBLE( A( J, L )*TEMP1+B( J, L )*TEMP2 )
  300.                   END IF
  301.   170          CONTINUE
  302.   180       CONTINUE
  303.          END IF
  304.       ELSE
  305. *
  306. *        Form  C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A +
  307. *                   C.
  308. *
  309.          IF( UPPER ) THEN
  310.             DO 210 J = 1, N
  311.                DO 200 I = 1, J
  312.                   TEMP1 = ZERO
  313.                   TEMP2 = ZERO
  314.                   DO 190 L = 1, K
  315.                      TEMP1 = TEMP1 + DCONJG( A( L, I ) )*B( L, J )
  316.                      TEMP2 = TEMP2 + DCONJG( B( L, I ) )*A( L, J )
  317.   190             CONTINUE
  318.                   IF( I.EQ.J ) THEN
  319.                      IF( BETA.EQ.DBLE( ZERO ) ) THEN
  320.                         C( J, J ) = DBLE( ALPHA*TEMP1+DCONJG( ALPHA )*
  321.      $                              TEMP2 )
  322.                      ELSE
  323.                         C( J, J ) = BETA*DBLE( C( J, J ) ) +
  324.      $                              DBLE( ALPHA*TEMP1+DCONJG( ALPHA )*
  325.      $                              TEMP2 )
  326.                      END IF
  327.                   ELSE
  328.                      IF( BETA.EQ.DBLE( ZERO ) ) THEN
  329.                         C( I, J ) = ALPHA*TEMP1 + DCONJG( ALPHA )*TEMP2
  330.                      ELSE
  331.                         C( I, J ) = BETA*C( I, J ) + ALPHA*TEMP1 +
  332.      $                              DCONJG( ALPHA )*TEMP2
  333.                      END IF
  334.                   END IF
  335.   200          CONTINUE
  336.   210       CONTINUE
  337.          ELSE
  338.             DO 240 J = 1, N
  339.                DO 230 I = J, N
  340.                   TEMP1 = ZERO
  341.                   TEMP2 = ZERO
  342.                   DO 220 L = 1, K
  343.                      TEMP1 = TEMP1 + DCONJG( A( L, I ) )*B( L, J )
  344.                      TEMP2 = TEMP2 + DCONJG( B( L, I ) )*A( L, J )
  345.   220             CONTINUE
  346.                   IF( I.EQ.J ) THEN
  347.                      IF( BETA.EQ.DBLE( ZERO ) ) THEN
  348.                         C( J, J ) = DBLE( ALPHA*TEMP1+DCONJG( ALPHA )*
  349.      $                              TEMP2 )
  350.                      ELSE
  351.                         C( J, J ) = BETA*DBLE( C( J, J ) ) +
  352.      $                              DBLE( ALPHA*TEMP1+DCONJG( ALPHA )*
  353.      $                              TEMP2 )
  354.                      END IF
  355.                   ELSE
  356.                      IF( BETA.EQ.DBLE( ZERO ) ) THEN
  357.                         C( I, J ) = ALPHA*TEMP1 + DCONJG( ALPHA )*TEMP2
  358.                      ELSE
  359.                         C( I, J ) = BETA*C( I, J ) + ALPHA*TEMP1 +
  360.      $                              DCONJG( ALPHA )*TEMP2
  361.                      END IF
  362.                   END IF
  363.   230          CONTINUE
  364.   240       CONTINUE
  365.          END IF
  366.       END IF
  367. *
  368.       RETURN
  369. *
  370. *     End of ZHER2K.
  371. *
  372.       END
  373.